Imports¶

In [ ]:
import pandas as pd
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats as stats
from itables import init_notebook_mode
import plotly.express as px


init_notebook_mode(all_interactive=True)

Loading data¶

In [ ]:
data = pd.read_excel(
    "../data/who_aap_2021_v9_11august2022.xlsx", sheet_name="AAP_2022_city_v9"
)
data.columns = [
    "WHO Region",
    "ISO3",
    "WHO Country Name",
    "City or Locality",
    "Measurement Year",
    "PM2.5",
    "PM10",
    "NO2",
    "PM2.5 temporal coverage",
    "PM10 temporal coverage",
    "NO2 temporal coverage",
    "Reference",
    "Number and type of monitoring stations",
    "Version of the database",
    "Status",
]
data
Out[ ]:
WHO Region ISO3 WHO Country Name City or Locality Measurement Year PM2.5 PM10 NO2 PM2.5 temporal coverage PM10 temporal coverage NO2 temporal coverage Reference Number and type of monitoring stations Version of the database Status
Loading... (need help?)

Biblioteka ITables 2.0 ułatwia wyszukiwanie konkretnej wartości w wyświetlanych danych (DataFrame'ach). Dodatkowo pozwala obejrzeć cały output, niezależnie od jego wielkości, co jest niewątpliwym plusem. W kontekście badanych danych używano jej np. do podglądania, czy dla danej wartości zaszła przewidywana zmiana, co się kryje w innej kolumnie wiersza, gdzie znajduje się szukana wartość oraz do wspomnianego wcześniej przeglądania dłuższych outputów.

Data review¶

Data types¶

In [ ]:
data.dtypes
Out[ ]:
0
Loading... (need help?)

Wszystkie typy danych są zgodne z logiką i prezentowaną zawartością.

Basic statistics¶

In [ ]:
data.describe()
Out[ ]:
Measurement Year PM2.5 PM10 NO2 PM2.5 temporal coverage PM10 temporal coverage NO2 temporal coverage Version of the database Status
Loading... (need help?)

Minimalne i maksymalne wartości w kolumnach "PM2.5 temporal coverage", "PM10 temporal coverage" oraz "NO2 temporal coverage" są zgodne z logiką (wiedząc, że są to procenty rocznego pokrycia pomiarami od 0 do 100), zatem nie ma w tych kolumnach wartości odstających. Kolumna "Measurement Year" i "Version of the database columns", bazując na przedstawionych podstawowych statystykach, również wydają się nie zawierać wartości odstających. Co do kolumn "PM2.5", "PM10" oraz "NO2" - mogą one zawierać wartości odstające, ale mając na uwadze, że są to dane dotyczące wartości zanieczyszczeń powietrza, nie zaleca się usuwania ich - zdecydowanie zaburzyłoby to dalszą analizę i ograniczyło jej sens.

NAs¶

In [ ]:
data.isna().sum()
Out[ ]:
0
Loading... (need help?)

Liczba wartości odstających (1) w kolumnie "WHO Region" budzi niepokój i zostanie to później sprawdzone. Wartości odstające w kolumnach "PM2.5", "PM10" i "NO2" wynikają z braku dokonania pomiaru w danym mieście w danym roku i ich typ można określić jako MAR - wartości tych brakuje często dla konkretnych lokalizacji lub lat. Braki w "PM2.5 temporal coverage", "PM10 temporal coverage" i "NO2 temporal coverage" również wynikają z braku wprowadzenia ich do bazy przez dane źródło, choć ich typ przypomina tu bardziej MCAR - nie zależy to od pozostałych danych. 5 NaNs w kolumnie "Reference" również zostanie zbadane w dalszej analizie. Kolumna "Number and type of monitoring stations" zawiera bardzo wiele wartości brakujących, ich typ również najbardziej przypomina MCAR. Kolumna "Status" w 100% składa się z NaNów, więc zostanie ona najprawdopodobniej usunięta na etapie selekcji atrybutów.

Checking the NA in WHO Region column¶

In [ ]:
data[data["WHO Region"].isna()]
Out[ ]:
WHO Region ISO3 WHO Country Name City or Locality Measurement Year PM2.5 PM10 NO2 PM2.5 temporal coverage PM10 temporal coverage NO2 temporal coverage Reference Number and type of monitoring stations Version of the database Status
Loading... (need help?)
In [ ]:
data["WHO Region"].unique()
Out[ ]:
array(['Eastern Mediterranean Region', 'European Region',
       'Region of the Americas', 'Western Pacific Region',
       'South East Asia Region', 'African Region', nan], dtype=object)
In [ ]:
na_slice = data[data["WHO Region"].isna()].copy()
na_slice.fillna({"WHO Region": "European Region"}, inplace=True)
data.update(na_slice)
data[data["ISO3"] == "LIE"]
Out[ ]:
WHO Region ISO3 WHO Country Name City or Locality Measurement Year PM2.5 PM10 NO2 PM2.5 temporal coverage PM10 temporal coverage NO2 temporal coverage Reference Number and type of monitoring stations Version of the database Status
Loading... (need help?)

Jedna brakująca wartość w kolumnie "WHO Region" istniała dla Lichtensteinu. Wynikało to najprawdopodobniej z błędu przy wprowadzaniu danych do bazy. Jako że Lichtenstein leży na terenie Europy założono, że należy on do tego samego regionu WHO co inne europejskie państwa, więc wypełniono brak odpowiednim regionem.

Checking NAs in Reference column¶

In [ ]:
data[data["Reference"].isna()]
Out[ ]:
WHO Region ISO3 WHO Country Name City or Locality Measurement Year PM2.5 PM10 NO2 PM2.5 temporal coverage PM10 temporal coverage NO2 temporal coverage Reference Number and type of monitoring stations Version of the database Status
Loading... (need help?)

Braki w kolumnie "Reference" występują tylko dla Dohy w Katarze. Z powodu niemożności oszacowania możliwego źródła danych z tego obszaru oraz mając na uwadze małe znaczenie tej informacji w kontekście przeprowadzanej analizy, zdecydowano pozostawić ten brak takim, jakim jest.

NAs in Number and type of monitoring stations column¶

Co do braków w kolumnie "Number and types of monitoring stations" liczba stacji może się zmieniać z czasem (dodawane są nowe, bądź aktualne mają awarię), więc uzupełnianie ich tą samą wartością, która znajduje się dla danego roku i miasta może wprowadzać w błąd. Dodatkowo, niektóre źródła nie zapewniają wcale takich informacji, więc biorąc to wszystko pod uwagę nic nie zrobiono w tymi brakami, jednocześnie wiedząc, że usunięcie ich znacznie skróciłoby dostępny zbiór danych. Typ tych braków można zakwalifikować jako MNAR, ponieważ może on zależeć od innych czynników niż dostępne.

Relation between variables¶

In [ ]:
sns.pairplot(data, hue="WHO Region").fig.suptitle(
    "Relations and distributions of variables vs WHO region", y=1.05, fontsize=20
)
plt.show()

Z przedstawionej figury można odczytać pewne relacje, np. dodatnia zależność liniowa między wartościami PM2.5 i PM10, co zostanie dogłębniej zbadane w dalszych etapach.

Feature selection and feature engineering¶

In [ ]:
data.drop(columns=["Version of the database", "Status"], inplace=True)

Zadecydowano o usunięciu kolumny "Status", ponieważ składała się ona w 100% w NaNów, więc nie wnosiła nic do analizy. Dodatkowo usunięto kolumnę "Version of the database", ponieważ nie wpływała ona na wyjaśnienie wartości zanieczyszczeń powietrza.

In [ ]:
def extract_and_sum_numbers(text):
    if pd.notna(text):
        numbers = [
            int(word.split()[0])
            for word in text.split(",")
            if word.strip()[0].isdigit()
        ]
        return sum(numbers)
    else:
        return np.nan


data["Total number of stations"] = data["Number and type of monitoring stations"].apply(
    extract_and_sum_numbers
)
data
Out[ ]:
WHO Region ISO3 WHO Country Name City or Locality Measurement Year PM2.5 PM10 NO2 PM2.5 temporal coverage PM10 temporal coverage NO2 temporal coverage Reference Number and type of monitoring stations Total number of stations
Loading... (need help?)

Dodano nową zmienną "Total number w stations", aby zbadać zależność ilości stacji pomiarowych od wartości zanieczyszczeń - czy np. bardziej zanieczyszczone regiony stawiają więcej stacji, aby być świadomym zagrożenia, czy też nie.

In [ ]:
region_dummies = pd.get_dummies(data["WHO Region"])
data_dummies = pd.concat([data, region_dummies], axis=1)
data_dummies
Out[ ]:
WHO Region ISO3 WHO Country Name City or Locality Measurement Year PM2.5 PM10 NO2 PM2.5 temporal coverage PM10 temporal coverage NO2 temporal coverage Reference Number and type of monitoring stations Total number of stations African Region Eastern Mediterranean Region European Region Region of the Americas South East Asia Region Western Pacific Region
Loading... (need help?)
In [ ]:
country_dummies = pd.get_dummies(data["WHO Country Name"])
data_dummies = pd.concat([data_dummies, country_dummies], axis=1)
data_dummies
Out[ ]:
WHO Region ISO3 WHO Country Name City or Locality Measurement Year PM2.5 PM10 NO2 PM2.5 temporal coverage PM10 temporal coverage NO2 temporal coverage Reference Number and type of monitoring stations Total number of stations African Region Eastern Mediterranean Region European Region Region of the Americas South East Asia Region Western Pacific Region Afghanistan Albania Algeria Andorra Argentina Australia Austria Bahamas Bahrain Bangladesh Belgium Bhutan Bolivia (Plurinational State of) Bosnia and Herzegovina Brazil Bulgaria Cameroon Canada Chile China Colombia Costa Rica Croatia Cuba Cyprus Czechia Denmark Ecuador Egypt El Salvador Estonia Ethiopia Fiji Finland France Georgia Germany Ghana Greece Guatemala Honduras Hungary Iceland India Indonesia Iran (Islamic Republic of) Iraq Ireland Israel Italy Jamaica Japan Jordan Kazakhstan Kenya Kuwait Kyrgyzstan Lao People's Democratic Republic Latvia Lebanon Liechtenstein Lithuania Luxembourg Madagascar Malaysia Maldives Malta Mauritius Mexico Monaco Mongolia Montenegro Morocco Myanmar Nepal Netherlands New Zealand Nigeria North Macedonia Norway Pakistan Panama Paraguay Peru Philippines Poland Portugal Qatar Republic of Korea Romania Russian Federation Saudi Arabia Senegal Serbia Singapore Slovakia Slovenia South Africa Spain Sri Lanka Sweden Switzerland Tajikistan Thailand Trinidad and Tobago Tunisia Turkey Turkmenistan Uganda Ukraine United Arab Emirates United Kingdom United Republic of Tanzania United States of America Uruguay Uzbekistan Venezuela (Bolivarian Republic of) Viet Nam
Loading... (need help?)

Do odpowiednich kopii ogólnego datasetu dodano także dummy values utworzone na podstawie kolumn "WHO Region" oraz "WHO Country Name", aby móc sprawdzić zależność między danym krajem lub regionem (do regresji punktowo-dwuseryjnej).

WHO Region-based analysis¶

In [ ]:
xlabels = [x.replace(" ", "\n") for x in data["WHO Region"].unique()]
xlabels

fig, ax = plt.subplots(2, 2, figsize=(18, 15))

sns.histplot(data=data, x="WHO Region", ax=ax[0, 0], shrink=0.9)
ax[0, 0].set_xticks(ticks=range(len(xlabels)), labels=xlabels)
ax[0, 0].set(title="Number of measurements vs WHO Region")

sns.boxplot(data=data, y="NO2", x="WHO Region", ax=ax[0, 1])
ax[0, 1].axhline(y=10, color="r")
ax[0, 1].set_xticks(ticks=range(len(xlabels)), labels=xlabels)
ax[0, 1].set(title="Distribution of NO2 values based on WHO Region for all years")

sns.boxplot(data=data, y="PM2.5", x="WHO Region", ax=ax[1, 0])
ax[1, 0].axhline(y=5, color="r")
ax[1, 0].set_xticks(ticks=range(len(xlabels)), labels=xlabels)
ax[1, 0].set(title="Distribution of PM2.5 values based on WHO Region for all years")

sns.boxplot(data=data, y="PM10", x="WHO Region", ax=ax[1, 1])
ax[1, 1].axhline(y=15, color="r")
ax[1, 1].set_xticks(ticks=range(len(xlabels)), labels=xlabels)
ax[1, 1].set(title="Distribution of PM10 values based on WHO Region for all years")

plt.suptitle("Data distribution based on WHO Region", fontsize=16)
plt.subplots_adjust(hspace=0.3, wspace=0.2)
plt.show()

Analizę rozpoczęto od generalizacji na regiony WHO. Wykres w lewym górnym rogu przedstawia ilość pomiarów wykonaną dla każdego regionu. Jak widać, zdecydowanie przoduje region europejski, co może świadczyć o największej świadomości ekologicznej w tym regionie. Najgorzej wypadają tu region Wschodniego Morza Śródziemnego oraz Afryka, głównie ze względu na to, że w tych regionach jest wiele krajów biednych, rozwijających się, do których świadomość ekologiczna jeszcze nie dotarła, bądź jest traktowana po macoszemu.\ Pozostałe wykresy przedstawiają rozkład wartości zanieczyszczeń z podziałem na regiony. Czerwona linia oznacza normę uznaną przez WHO dla każdego z zanieczyszczeń. Pokazuje to jak złym powietrzem oddycha się na całym świecie. Najlepiej pod względem PM2.5 wypada region Ameryk oraz Europa, mając jednak na uwadze liczne wartości powyżej 95 percentyla. Najgorsza sytuacja panuje w regionie Zachodniego Pacyfiku (gdzie normy zawyżają Chiny oraz Wietnam) oraz w Azji Południowo-Wschodniej (tutaj głównym winowajcą są Indie). W przypadku PM10 sytuacja wygląda podobnie pod względem niskich wartości - tu przodują wciąż Europa i Ameryki, ale tym razem najgorsze powietrze (pod względem PM10) jest w regionie Wschodniego Morza Śródziemniego oraz w Azji Południowo-Wschodniej. Analogicznie badając wartości NO2.

Correlation¶

In [ ]:
data.columns
Out[ ]:
Index(['WHO Region', 'ISO3', 'WHO Country Name', 'City or Locality',
       'Measurement Year', 'PM2.5', 'PM10', 'NO2', 'PM2.5 temporal coverage',
       'PM10 temporal coverage', 'NO2 temporal coverage', 'Reference',
       'Number and type of monitoring stations', 'Total number of stations'],
      dtype='object')
In [ ]:
corr_cols = [
    "Measurement Year",
    "PM2.5",
    "PM10",
    "NO2",
    "PM2.5 temporal coverage",
    "PM10 temporal coverage",
    "NO2 temporal coverage",
    "Total number of stations",
]

corr_df = data[corr_cols].corr()

Wybrano jedynie takie kolumny do analizy korelacji liniowej Pearsona, które zawierają wartości ciągłe.

In [ ]:
fig, ax = plt.subplots(1, 1, figsize=(6, 6))

sns.heatmap(
    corr_df,
    ax=ax,
    vmin=-1,
    vmax=1,
    center=0,
    annot=True,
    fmt=".2f",
    cmap="viridis",
    annot_kws={"fontsize": 10},
)
ax.set(title="Correlation heatmap for the whole dataset")

plt.show()

Wartości korelacji liniowej dla całego zbioru danych wyraźnie wskazują istotną zależność między wartościami PM2.5 oraz PM10, a także mniejsze pomiędzy NO2 a PM2.5 oraz PM10. Ponadto można odczytać niewielkie ujemne zależności między wartościami PM2.5, PM10 i NO2 a stopniem rocznego pokrycia każdej z nich - jest to ciekawe, ponieważ pokazuje, że im częściej wykonywano pomiary, tym mniejsze były zanieczyszczenia, co może być pośrednio z korelowane z tym, że większe pokrycie występowało najprawdopodobniej w krajach bardziej rozwiniętych. Ponadto wystepują dodatnie zależności między pokryciami każdego ze współczynników, co jest dość logiczne - jeśli stacja mierzyła jeden, to najczęściej i pozostałe współczynniki, i na odwrót. Występuje też niewielka dodatnia zależność pomiędzy ilością stacji pomiarowych a wartością NO2.

In [ ]:
i = 0
j = 0

fig, ax = plt.subplots(2, 3, figsize=(23, 15))

for r in data["WHO Region"].unique():
    corr_df = data[data["WHO Region"] == r][corr_cols].corr()
    sns.heatmap(
        corr_df,
        ax=ax[i][j],
        vmin=-1,
        vmax=1,
        center=0,
        annot=True,
        fmt=".2f",
        cmap="viridis",
        annot_kws={"fontsize": 10},
    )
    ax[i][j].set(title=f"Correlation heatmap for the {r}")
    i += 1
    i %= 2
    j += 1
    j %= 3

plt.subplots_adjust(hspace=0.5, wspace=0.5)
plt.suptitle("Correlation heatmaps for an each region", fontsize=16)
plt.show()

Na powyższych heatmapach przedstawiono zależności pomiędzy poszczególnymi zmiennymi dla każdego z regionów. Wartości współczynników korelacji nieraz mocno różnią się od siebie, jednak wszędzie występują duże dodatnie zależności pomiędzy wartościami każdego ze współczynników, szczególnie PM2.5 i PM10. Ponadto w regionie Azji Południowo-Wschodniej występuje ujemna zależność między rokiem pomiaru a wartością PM2.5. Wskazuje to na nieznaczną poprawę jakości powietrza w tym regionie wraz z upływem lat. Dla regionu Ameryk wyróżniają się ujemne zależności między wartościami współczynników a ich stopniem rocznego pokrycia. W regionie europejskim wyróżniają się wysokie wartości współczynników dla stopni pokrycia dla siebie nawzajem, co może wskazywać, że zazwyczaj europejskie czujniki mierzą wszystkie z 3 rodzajów zanieczyszczeń. Dla Afryki widać, że zupełnie nie wypełniano kolumny z rocznym pokryciem pomiarami współczynnika NO2 oraz, że występuje bardzo silna pozytywna zależność między stopniem pokrycia PM2.5 i PM10.

Point-biserial correlation¶

In [ ]:
for r in data_dummies["WHO Region"].unique():
    print(f"\n{r}:")
    for c in corr_cols:
        temp_df = data_dummies.dropna(subset=c)
        print(f"{r} vs {c}: {stats.pointbiserialr(temp_df[r], temp_df[c])}")
Eastern Mediterranean Region:
Eastern Mediterranean Region vs Measurement Year: SignificanceResult(statistic=-0.020433493981191153, pvalue=0.00024600122593072257)
Eastern Mediterranean Region vs PM2.5: SignificanceResult(statistic=0.1031366363946883, pvalue=7.199182966700921e-37)
Eastern Mediterranean Region vs PM10: SignificanceResult(statistic=0.36426022043226786, pvalue=0.0)
Eastern Mediterranean Region vs NO2: SignificanceResult(statistic=0.17002278268537616, pvalue=1.264128529846678e-143)
Eastern Mediterranean Region vs PM2.5 temporal coverage: SignificanceResult(statistic=-0.14136809403182024, pvalue=8.62450943600881e-34)
Eastern Mediterranean Region vs PM10 temporal coverage: SignificanceResult(statistic=-0.07623380154639861, pvalue=2.1558817653567293e-08)
Eastern Mediterranean Region vs NO2 temporal coverage: SignificanceResult(statistic=-0.007052246771967366, pvalue=0.31995921477155914)
Eastern Mediterranean Region vs Total number of stations: SignificanceResult(statistic=0.0350995745864689, pvalue=0.0010187075804140191)

European Region:
European Region vs Measurement Year: SignificanceResult(statistic=-0.06181305313485798, pvalue=1.2461816552351145e-28)
European Region vs PM2.5: SignificanceResult(statistic=-0.5019660662896209, pvalue=0.0)
European Region vs PM10: SignificanceResult(statistic=-0.4709508188792565, pvalue=0.0)
European Region vs NO2: SignificanceResult(statistic=-0.034143428574115814, pvalue=3.609955837755273e-07)
European Region vs PM2.5 temporal coverage: SignificanceResult(statistic=0.08638815294114853, pvalue=1.5723213252180334e-13)
European Region vs PM10 temporal coverage: SignificanceResult(statistic=0.3171493900311768, pvalue=5.152144759387591e-126)
European Region vs NO2 temporal coverage: SignificanceResult(statistic=0.36944088974020795, pvalue=0.0)
European Region vs Total number of stations: SignificanceResult(statistic=-0.02608689386184582, pvalue=0.01463093702593799)

Region of the Americas:
Region of the Americas vs Measurement Year: SignificanceResult(statistic=-0.08773089342540233, pvalue=4.977480469986595e-56)
Region of the Americas vs PM2.5: SignificanceResult(statistic=-0.26049484419636904, pvalue=6.489615190371197e-232)
Region of the Americas vs PM10: SignificanceResult(statistic=0.0017056317906980131, pvalue=0.8042925071377917)
Region of the Americas vs NO2: SignificanceResult(statistic=-0.13093600997685742, pvalue=1.7978730079789976e-85)
Region of the Americas vs PM2.5 temporal coverage: SignificanceResult(statistic=0.05833064553479305, pvalue=6.404069137069385e-07)
Region of the Americas vs PM10 temporal coverage: SignificanceResult(statistic=-0.08615753766238893, pvalue=2.4449908290095345e-10)
Region of the Americas vs NO2 temporal coverage: SignificanceResult(statistic=-0.000721393180807116, pvalue=0.9189685457267649)
Region of the Americas vs Total number of stations: SignificanceResult(statistic=-0.10006409915353032, pvalue=6.1940915210474116e-21)

Western Pacific Region:
Western Pacific Region vs Measurement Year: SignificanceResult(statistic=0.1907110651035102, pvalue=2.5510330156250427e-261)
Western Pacific Region vs PM2.5: SignificanceResult(statistic=0.6279332135664273, pvalue=0.0)
Western Pacific Region vs PM10: SignificanceResult(statistic=0.01417317443144923, pvalue=0.039475752565946924)
Western Pacific Region vs NO2: SignificanceResult(statistic=0.15804727736142404, pvalue=3.940644188511663e-124)
Western Pacific Region vs PM2.5 temporal coverage: SignificanceResult(statistic=0.1547084989920406, pvalue=3.3040342092080696e-40)
Western Pacific Region vs PM10 temporal coverage: SignificanceResult(statistic=0.058758057555798396, pvalue=1.6109785653262883e-05)
Western Pacific Region vs NO2 temporal coverage: SignificanceResult(statistic=-0.1327902656663876, pvalue=6.272854169555007e-79)
Western Pacific Region vs Total number of stations: SignificanceResult(statistic=0.026340006579452183, pvalue=0.013698117063235612)

South East Asia Region:
South East Asia Region vs Measurement Year: SignificanceResult(statistic=-0.030471870781641, pvalue=4.5443664453860183e-08)
South East Asia Region vs PM2.5: SignificanceResult(statistic=0.21706787449153575, pvalue=6.106384672875838e-160)
South East Asia Region vs PM10: SignificanceResult(statistic=0.5869955522796655, pvalue=0.0)
South East Asia Region vs NO2: SignificanceResult(statistic=0.040576068267698405, pvalue=1.467932417415133e-09)
South East Asia Region vs PM2.5 temporal coverage: SignificanceResult(statistic=-0.34018669806329116, pvalue=1.585732162934143e-196)
South East Asia Region vs PM10 temporal coverage: SignificanceResult(statistic=-0.25030145383864205, pvalue=1.1634152500678984e-77)
South East Asia Region vs NO2 temporal coverage: SignificanceResult(statistic=-0.3539678931236533, pvalue=0.0)
South East Asia Region vs Total number of stations: SignificanceResult(statistic=0.08215839173536814, pvalue=1.3537639103873756e-14)

African Region:
African Region vs Measurement Year: SignificanceResult(statistic=0.0165086116259789, pvalue=0.0030561082394488653)
African Region vs PM2.5: SignificanceResult(statistic=0.02214265459012913, pvalue=0.006600545585984725)
African Region vs PM10: SignificanceResult(statistic=0.0637242950419178, pvalue=1.908574190231089e-20)
African Region vs NO2: SignificanceResult(statistic=0.015695202706818095, pvalue=0.019359000014242633)
African Region vs PM2.5 temporal coverage: SignificanceResult(statistic=-0.16552942878503737, pvalue=7.465659578923453e-46)
African Region vs PM10 temporal coverage: SignificanceResult(statistic=-0.19710229395003298, pvalue=2.8680642200596e-48)
African Region vs NO2 temporal coverage: SignificanceResult(statistic=-0.1511563002336828, pvalue=5.683997130702588e-102)
African Region vs Total number of stations: SignificanceResult(statistic=0.027467689198414545, pvalue=0.010150509323086856)

Badając wartości współczynników korelacji punktowo-dwuseryjnej zaobserowowano:

  • dla regionu Wschodniego Morza Śródziemnego występuje dodatnia zależność wobec wartości PM10 - czyli tutaj są one zazwyczaj większe niż w innych regionach
  • w Europie występują wysokie ujemne zależności wobec PM2.5 oraz PM10 - powietrze jest nieco czystsze
  • dla regionu Ameryk występuje niewielka ujemna zależność wobec PM2.5
  • dla Azji Południowo-Wschodniej zachodzi duża dodatnia zależność wobec wartości PM10 - w tym regionie są one wyższe niż w pozostałych; ponadto występuje tu zazwyczaj wyższe roczne pokrycie pomiarami
  • dla Afryki nie zauważono żadnych znaczących zależności

Air pollution across years¶

In [ ]:
grouped_df_pm25 = data.groupby(["WHO Region", "Measurement Year"])["PM2.5"].mean(
    numeric_only=True
)
grouped_df_pm10 = data.groupby(["WHO Region", "Measurement Year"])["PM10"].mean(
    numeric_only=True
)
grouped_df_no2 = data.groupby(["WHO Region", "Measurement Year"])["NO2"].mean(
    numeric_only=True
)
In [ ]:
fig, ax = plt.subplots(3, 1, figsize=(15, 12))

sns.lineplot(grouped_df_pm25, x="Measurement Year", y=grouped_df_pm25.values, ax=ax[0])
ax[0].axhline(y=5, color="r")
ax[0].set(title="PM2.5 pullution", ylabel="PM2.5 [(μg/m3]")

sns.lineplot(grouped_df_pm10, x="Measurement Year", y=grouped_df_pm10.values, ax=ax[1])
ax[1].axhline(y=15, color="r")
ax[1].set(title="PM10 pullution", ylabel="PM10 [(μg/m3]")

sns.lineplot(grouped_df_no2, x="Measurement Year", y=grouped_df_no2.values, ax=ax[2])
ax[2].axhline(y=10, color="r")
ax[2].set(title="NO2 pullution", ylabel="NO2 [(μg/m3]")

plt.suptitle("Air pollution across years for all data", fontsize=16)
plt.subplots_adjust(hspace=0.3)
plt.show()

Dla całości danych przedstawiono wartości współczynników zanieczyszczenia powietrza dla kolejnych lat. Jasnoniebieska obwoluta pokazuje przedziały ufności, a raczej w tym przypadku - maksymalne i minimalne wartości średnich dla każdego z regionów. Czerwoną linią oznaczona normę przyjętą przez WHO dla każdego ze współczynników. W przypadku każdego z nich widać generalną (delikatną) tendencję spadkową w ostatnich latach, jednak z przyrostem w ostatnim badanym roku dla PM2.5 oraz PM10.

In [ ]:
i = 0
j = 0

fig, ax = plt.subplots(3, 2, figsize=(15, 12))

for r in data["WHO Region"].unique():

    sub_grouped_df = grouped_df_pm25[
        grouped_df_pm25.index.get_level_values("WHO Region").isin([r])
    ]

    sns.lineplot(
        sub_grouped_df,
        x="Measurement Year",
        y=sub_grouped_df.values,
        ax=ax[i, j],
    )
    ax[i, j].axhline(y=5, color="r")
    ax[i, j].set(title=f"PM2.5 pullution for {r}", ylabel="PM2.5 [(μg/m3]")
    i += 1
    i %= 3
    j += 1
    j %= 2

plt.suptitle("PM2.5 pollution across years for all regions", fontsize=16)
plt.subplots_adjust(hspace=0.3)
plt.show()

Badając samo PM2.5 dla każdego z regionów można wykryć wyraźny ujemny trend w regionach Azji Południowo-Wschodniej, Europie, Amerykach, a w ostatnich latach także w regionie Zachodniego Pacyfiku. W regionie afrykańskim wartości zmieniają się w sposób dość losowy. Niestety, w każdym z regionów norma jest znacząco przekroczona (oprócz ostatniego roku w Amerykach). Zazwyczaj, średnie wartości zanieczyszczeń PM2.5 znacznie różniły się między regionami na korzyść Europy i Ameryk.

In [ ]:
i = 0
j = 0

fig, ax = plt.subplots(3, 2, figsize=(15, 12))

for r in data["WHO Region"].unique():

    sub_grouped_df = grouped_df_pm10[
        grouped_df_pm10.index.get_level_values("WHO Region").isin([r])
    ]

    sns.lineplot(
        sub_grouped_df,
        x="Measurement Year",
        y=sub_grouped_df.values,
        ax=ax[i, j],
    )
    ax[i, j].axhline(y=15, color="r")
    ax[i, j].set(title=f"PM10 pullution for {r}", ylabel="PM10 [(μg/m3]")
    i += 1
    i %= 3
    j += 1
    j %= 2

plt.suptitle("PM10 pollution across years for all regions", fontsize=16)
plt.subplots_adjust(hspace=0.3)
plt.show()

W przypadku PM10 ujemny trend występuje jedynie w Europie oraz Amerykach. Poza tymi regionami utrzymuje się on na mniej więcej stałym poziomie, ze znaczącym skokiem w okolicach roku 2009 dla regionu Wchodniego Morza Śródziemnego. Norma nie została średnio przekroczona tylko w ostatnich latach dla Europy i Ameryk. Warto zwrócić uwagę na podziałkę osi Y. W rejonach Wschodniego Morza Śródziemnego oraz Azji Południowo-Wschodniej zanieczyszczenie przyjmowało ogrmne wartości w porównaniu z pozostałymi.

In [ ]:
i = 0
j = 0

fig, ax = plt.subplots(3, 2, figsize=(15, 12))

for r in data["WHO Region"].unique():

    sub_grouped_df = grouped_df_no2[
        grouped_df_no2.index.get_level_values("WHO Region").isin([r])
    ]

    sns.lineplot(
        sub_grouped_df,
        x="Measurement Year",
        y=sub_grouped_df.values,
        ax=ax[i, j],
    )
    ax[i, j].axhline(y=10, color="r")
    ax[i, j].set(title=f"NO2 pullution for {r}", ylabel="NO2 [(μg/m3]")
    i += 1
    i %= 3
    j += 1
    j %= 2

plt.suptitle("NO2 pollution across years for all regions", fontsize=16)
plt.subplots_adjust(hspace=0.3)
plt.show()

Wartości NO2 charakteryzuje niewielki ujemny trend dla każdego z regionów oprócz Azji Południowo-Wschodniej, gdzie istnieje niewielki trend wzrostowy. Najwyższe średnie wartości były rejestrowane w rejonie Wschodniego Morza Śródziemnego.

Most and least polluted regions¶

PM2.5¶
In [ ]:
grouped_region_pm25 = data.groupby("WHO Region")["PM2.5"].mean()
pd.DataFrame(grouped_region_pm25).sort_values(by="PM2.5")
Out[ ]:
PM2.5
WHO Region
Loading... (need help?)

Zgodnie z poprzednimi spostrzeżeniami, najniższe średnie zanieczyszczenie PM2.5 występuje w regionie Ameryk oraz Europy, a najwyższe w rejonie Zachodniego Pacyfiku oraz Azji Południowo-Wschodniej.

In [ ]:
grouped_region_pm25 = data.groupby("WHO Region")["PM2.5"].median()
pd.DataFrame(grouped_region_pm25).sort_values(by="PM2.5")
Out[ ]:
PM2.5
WHO Region
Loading... (need help?)

Mediana zanieczyszczeń PM2.5 nie odbiega znacząco od średniej, co oznacza, że całe regiony są równomierne zanieczyszczone PM2.5.

PM10¶
In [ ]:
grouped_region_pm10 = data.groupby("WHO Region")["PM10"].mean()
pd.DataFrame(grouped_region_pm10).sort_values(by="PM10")
Out[ ]:
PM10
WHO Region
Loading... (need help?)

Rejony najmniej zanieczyszczone cząsteczkami PM10 są analogiczne do tych wspomianych w przypadku PM2.5. Jednakże tutaj Europa ma pod tym względem lepsze powietrze niż Ameryki, które średnio tylko odrobinę przegrywają z Zachodnim Pacyfikiem, który miał najgorszy wynik PM2.5. Zdecydowanie najbardziej zanieczyszczonym regionem PM10 jest Wschodnie Morze Śródziemne.

In [ ]:
grouped_region_pm10 = data.groupby("WHO Region")["PM10"].median()
pd.DataFrame(grouped_region_pm10).sort_values(by="PM10")
Out[ ]:
PM10
WHO Region
Loading... (need help?)

Wyniki bazujące na medianie nieco odbiegają od tych opartych o średnią, co wskazują na istotne różnice w zanieczyszczeniach wewnątrz regionów z pojedynczymi krajami/miastami znacząco zanieczyszczonymi.

NO2¶
In [ ]:
grouped_region_no2 = data.groupby("WHO Region")["NO2"].mean()
pd.DataFrame(grouped_region_no2).sort_values(by="NO2")
Out[ ]:
NO2
WHO Region
Loading... (need help?)

Najczystsze regiony pod względem NO2 to również Europa i Ameryki, a najbardziej zanieczyszczona - Zachodni Pacyfik oraz Wschodnie Morze Śródziemne.

In [ ]:
grouped_region_no2 = data.groupby("WHO Region")["NO2"].median()
pd.DataFrame(grouped_region_no2).sort_values(by="NO2")
Out[ ]:
NO2
WHO Region
Loading... (need help?)

Mediana prezentuje się tutaj bardzo podobnie do średniej, tak jak w przypadku PM2.5.

Stations analysis¶

In [ ]:
grouped_region_no_stations = data.groupby("WHO Region")[
    "Total number of stations"
].mean()
pd.DataFrame(grouped_region_no_stations).sort_values(by="Total number of stations")
Out[ ]:
Total number of stations
WHO Region
Loading... (need help?)

Co ciekawe najmniejszą średnią liczbę stacji posiada region Ameryk, a największą - region Wschodniego Morza Śródziemnego. Należy mieć tu jednak na uwadze liczne braki w tej kolumnie.

In [ ]:
grouped_region_no_stations = data.groupby("WHO Region")[
    "Total number of stations"
].median()
pd.DataFrame(grouped_region_no_stations).sort_values(by="Total number of stations")
Out[ ]:
Total number of stations
WHO Region
Loading... (need help?)

Mediana wyników prezentuje się już bardziej jednostajnie, wciąż mniej więcej utrzymując kolejność regionów.

Country-based analysis¶

In [ ]:
fig, ax = plt.subplots(2, 2, figsize=(40, 32))

sns.histplot(data=data, x="WHO Country Name", ax=ax[0, 0], shrink=0.7)
ax[0, 0].set_xticklabels(ax[0, 0].get_xticklabels(), rotation=90)
ax[0, 0].set(title="Number of measurements vs Country")

sns.boxplot(data=data, y="NO2", x="WHO Country Name", ax=ax[0, 1], width=0.7)
ax[0, 1].axhline(y=10, color="r")
ax[0, 1].grid(axis="x")
ax[0, 1].set_xticklabels(ax[0, 1].get_xticklabels(), rotation=90)
ax[0, 1].set(title="Distribution of NO2 values based on Country for all years")

sns.boxplot(data=data, y="PM2.5", x="WHO Country Name", ax=ax[1, 0], width=0.7)
ax[1, 0].axhline(y=5, color="r")
ax[1, 0].grid(axis="x")
ax[1, 0].set_xticklabels(ax[1, 0].get_xticklabels(), rotation=90)
ax[1, 0].set(title="Distribution of PM2.5 values based on Country for all years")

sns.boxplot(data=data, y="PM10", x="WHO Country Name", ax=ax[1, 1], width=0.7)
ax[1, 1].axhline(y=15, color="r")
ax[1, 1].grid(axis="x")
ax[1, 1].set_xticklabels(ax[1, 1].get_xticklabels(), rotation=90)
ax[1, 1].set(title="Distribution of PM10 values based on Country for all years")

plt.suptitle("Data distribution based on Country", fontsize=20)
plt.subplots_adjust(hspace=0.3)
plt.show()
C:\Users\adamp\AppData\Local\Temp\ipykernel_40448\2476146839.py:4: UserWarning: FixedFormatter should only be used together with FixedLocator
  ax[0, 0].set_xticklabels(ax[0, 0].get_xticklabels(), rotation=90)

Przechodząc do analizy opartej o kraje, ponownie rozpoczęto od zbadania rozkładu danych. Z histogramu w lewym górnym rogu widać, że w liczbie pomiarów zdecydowanie przodują Chiny znane z ogromnego zanieczyszczenia powietrza. Zaraz za nimi: (co ciekawe) Włochy, Niemcy, Hiszpania i Indie.\ Pod względem rozkładu wartości zanieczyszczeń PM2.5 można wyróżnić kraje z niewielkimi IQR jak np. Grecja, Malta czy USA, a także z bardzo dużymi, towarzyszącej im przy tym wysokiej medianie, tutaj: Bangladesz, Chiny, Kamerun, Mongolia, Pakista, Senegal. Chiny oraz Indie mają przy tym bardzo wiele wartości odstających powyżej 95 percentyla, co sugeruje, że istnieje wiele miast z bardzo wysokim zanieczyszczeniem PM2.5.\ Badając rozkład PM10 uwagę zwraca na siebie szczególnie Pakistan z ogromnym IQR oraz bardzo wysoką medianą. Ponadto duże zanieczyszczenia mają także: Bahrajn, Egipt, Ghana, Irak, Mongolia i Katar.\ W przypadku NO2 wartości są już mniej zróżnicowane, jednak przoduje tutaj Iran, Irak oraz Bahrajn - wszystkie położone geograficznie bardzo blisko siebie.

Most and least polluted countries (across all years)¶

Pokuszono się o wyspecyfikowanie najmniej oraz najbardziej zanieczyszczonych państw pod względem każdego ze współczynników, zarówno pod wzgledem średniej oraz mediany.

PM2.5¶
In [ ]:
grouped_country_pm25 = data.groupby("WHO Country Name")["PM2.5"].mean()
grouped_country_pm25 = (
    pd.DataFrame(grouped_country_pm25).sort_values(by="PM2.5").dropna()
)
grouped_country_pm25.index[:3].to_list()
Out[ ]:
['Bahamas', 'Iceland', 'Estonia']
In [ ]:
grouped_country_pm25.index[-3:].to_list()
Out[ ]:
['Bangladesh', 'Cameroon', 'Afghanistan']
In [ ]:
grouped_country_pm25 = data.groupby("WHO Country Name")["PM2.5"].median()
grouped_country_pm25 = (
    pd.DataFrame(grouped_country_pm25).sort_values(by="PM2.5").dropna()
)
grouped_country_pm25.index[:3].to_list()
Out[ ]:
['Bahamas', 'Iceland', 'Estonia']
In [ ]:
grouped_country_pm25.index[-3:].to_list()
Out[ ]:
['Tajikistan', 'Mongolia', 'Afghanistan']

Mediana nie wprowadziła tu znacznych zmian w stosunku do średniej, przez co można uznać Bahamy, Islandię i Estonię za najmniej zanieczyszczone kraje cząsteczkami PM2.5, a Afganistan wraz z towarzyszącymi mu na listach państwami - za najbardziej.

PM10¶
In [ ]:
grouped_country_pm10 = data.groupby("WHO Country Name")["PM10"].mean()
grouped_country_pm10 = (
    pd.DataFrame(grouped_country_pm10).sort_values(by="PM10").dropna()
)
grouped_country_pm10.index[:3].to_list()
Out[ ]:
['Bahamas', 'Iceland', 'Estonia']
In [ ]:
grouped_country_pm10.index[-3:].to_list()
Out[ ]:
['Ghana', 'Egypt', 'Pakistan']
In [ ]:
grouped_country_pm10 = data.groupby("WHO Country Name")["PM10"].median()
grouped_country_pm10 = (
    pd.DataFrame(grouped_country_pm10).sort_values(by="PM10").dropna()
)
grouped_country_pm10.index[:3].to_list()
Out[ ]:
['Bahamas', 'Iceland', 'Finland']
In [ ]:
grouped_country_pm25.index[-3:].to_list()
Out[ ]:
['Tajikistan', 'Mongolia', 'Afghanistan']

Pod względem PM10 również najlepiej wypadają Bahamy, Islandia ora Estonia, a biorąc pod uwagę medianę, także Finlandia. Najbardziej zanieczyszczone kraje to ponownie: Afganistan, Mongolia, Tadżykistan oraz Ghana i Egipt.

NO2¶
In [ ]:
grouped_country_no2 = data.groupby("WHO Country Name")["NO2"].mean()
grouped_country_no2 = pd.DataFrame(grouped_country_no2).sort_values(by="NO2").dropna()
grouped_country_no2.index[:3].to_list()
Out[ ]:
['Estonia', 'Myanmar', 'Iceland']
In [ ]:
grouped_country_no2.index[-3:].to_list()
Out[ ]:
['Lebanon', 'Iraq', 'Iran (Islamic Republic of)']
In [ ]:
grouped_country_no2 = data.groupby("WHO Country Name")["NO2"].median()
grouped_country_no2 = pd.DataFrame(grouped_country_no2).sort_values(by="NO2").dropna()
grouped_country_no2.index[:3].to_list()
Out[ ]:
['New Zealand', 'Iceland', 'Estonia']
In [ ]:
grouped_country_no2.index[-3:].to_list()
Out[ ]:
['Lebanon', 'Iraq', 'Iran (Islamic Republic of)']

Badając wartości NO2, ponownie najczystsza jest Estonia i Islandia, a dodatkowo Birma i Nowa Zelandia, a najbardziej zanieczyszczone: Iran, Irak i Liban.

Biggest progress and regress in air polution¶

Zbadano, czy na przestrzeni lat jakieś państwa odnotowały postęp w zakresie czystego powietrza, czy też wręcz przeciwne. Dodatnia wartość różnicy oznacza progress, a ujemna - regress.

PM2.5¶
In [ ]:
grouped = data.groupby(["WHO Country Name", "Measurement Year"]).mean(numeric_only=True)
grouped_pm25 = grouped.dropna(subset="PM2.5")

earliest = grouped_pm25.groupby(level=0)["PM2.5"].first()
latest = grouped_pm25.groupby(level=0)["PM2.5"].last()

difference = earliest - latest  # + --> progress; - --> regress

difference
Out[ ]:
PM2.5
WHO Country Name
Loading... (need help?)
In [ ]:
difference.sort_values()[:3]
Out[ ]:
PM2.5
WHO Country Name
Loading... (need help?)
In [ ]:
difference.sort_values()[-3:]
Out[ ]:
PM2.5
WHO Country Name
Loading... (need help?)

Jak widać największy regres pod względem PM2.5 odnotowały: Indonezja, Madagaskar oraz Mongolia, a progres: Wietnam, Uganda i Pakistan.

PM10¶
In [ ]:
grouped_pm10 = grouped.dropna(subset="PM10")

earliest = grouped_pm10.groupby(level=0)["PM10"].first()
latest = grouped_pm10.groupby(level=0)["PM10"].last()

difference = earliest - latest  # + --> progress; - --> regres

difference
Out[ ]:
PM10
WHO Country Name
Loading... (need help?)
In [ ]:
difference.sort_values()[:3]
Out[ ]:
PM10
WHO Country Name
Loading... (need help?)
In [ ]:
difference.sort_values()[-3:]
Out[ ]:
PM10
WHO Country Name
Loading... (need help?)

Pod względem cząsteczek PM10 ogromny regres zarejestrowano w: Bhutanie, Egipcie i na Madagaskarze, a równie duży progres w: Chinach, Iraku i Pakistanie.

NO2¶
In [ ]:
grouped_no2 = grouped.dropna(subset="NO2")

earliest = grouped_no2.groupby(level=0)["NO2"].first()
latest = grouped_no2.groupby(level=0)["NO2"].last()

difference = earliest - latest  # + --> progress; - --> regres

difference
Out[ ]:
NO2
WHO Country Name
Loading... (need help?)
In [ ]:
difference.sort_values()[:3]
Out[ ]:
NO2
WHO Country Name
Loading... (need help?)
In [ ]:
difference.sort_values()[-3:]
Out[ ]:
NO2
WHO Country Name
Loading... (need help?)

W przypadku wartości NO2 regres odnotowano w: Meksyku, Iranie i Zjednoczonych Emiratach Arabskich, a postęp w: Litwie, Południowej Afryce oraz Kuwejcie.

Stations analysis¶

In [ ]:
summed_stations_df = data.groupby("WHO Country Name")["Total number of stations"].sum()
mean_stations_df = data.groupby("WHO Country Name")["Total number of stations"].mean()

fig, ax = plt.subplots(2, 1, figsize=(20, 15))

sns.barplot(
    x=summed_stations_df.index, y=summed_stations_df.values, ax=ax[0], color="gray"
)
ax[0].set_xticklabels(ax[0].get_xticklabels(), rotation=90)
ax[0].set(
    title="Total number of station per country", ylabel="Summed number of stations"
)

sns.barplot(
    x=mean_stations_df.index, y=mean_stations_df.values, ax=ax[1], color="darkgrey"
)
ax[1].set_xticklabels(ax[1].get_xticklabels(), rotation=90)
ax[1].set(title="Mean number of station per country", ylabel="Mean number of stations")

plt.suptitle("Stations distribution", fontsize=16)
plt.subplots_adjust(hspace=0.7)
plt.show()

Podjęto się również analizy ilości stacji pomiarowych w każdym z krajów, zarówno pod względem ich całkowitej ilości, jak i średniej. W sumarycznej ilości wyróżniają się Indie, a poza nimi, choć znacznie mniej: Szwajcaria, Włochy, Niemcy, Czechy, Kanada, Francja i Hiszpania. W średniej ilości stacji prowadzi za to Gwatemala, a zaraz za nią są: Iran, Zjednoczone Emiraty Arabskie oraz Urugwaj. Należy pamiętać, że dane dotyczące ilości stacji mogą mieć wysoki bias, ze względu na dużą ilość wartości brakujących typu MCAR.

Air pollution map¶

PM2.5¶
In [ ]:
grouped_country_pm25 = data.groupby(["WHO Country Name", "ISO3"])["PM2.5"].mean()
df_grouped_country_pm25 = grouped_country_pm25.reset_index()

fig = px.choropleth(
    df_grouped_country_pm25,
    locations="ISO3",
    locationmode="ISO-3",
    color="PM2.5",
    hover_name="WHO Country Name",
    color_continuous_scale=px.colors.sequential.Sunsetdark,
)

fig.update_layout(
    title_text="Average concentration of PM2.5 across all years",
)

fig.show()

Mapa świata z naniesionymi wartościami zanieczyszczeń PM2.5 wyraźnie przedstawia całą Azję Południową jako nabardziej zanieczyszczoną. Najlepiej prezentuje się Ameryka północna oraz Europa Północna i Zachodnia.

PM10¶
In [ ]:
grouped_country_pm10 = data.groupby(["WHO Country Name", "ISO3"])["PM10"].mean()
df_grouped_country_pm10 = grouped_country_pm10.reset_index()

fig = px.choropleth(
    df_grouped_country_pm10,
    locations="ISO3",
    locationmode="ISO-3",
    color="PM10",
    hover_name="WHO Country Name",
    color_continuous_scale=px.colors.sequential.Sunsetdark,
)

fig.update_layout(
    title_text="Average concentration of PM10 across all years",
)

fig.show()

Podobnie jak mapa dla PM2.5, mapa dla PM10 pokazuje te same siedliska zanieczyszczeń, dodając przy tym niektóre państwa afrykańskie. Należy zauważyć, jak przerażająco prezentuje się ona dla Pakistanu.

NO2¶
In [ ]:
grouped_country_no2 = data.groupby(["WHO Country Name", "ISO3"])["NO2"].mean()
df_grouped_country_no2 = grouped_country_no2.reset_index()

fig = px.choropleth(
    df_grouped_country_no2,
    locations="ISO3",
    locationmode="ISO-3",
    color="NO2",
    hover_name="WHO Country Name",
    color_continuous_scale=px.colors.sequential.Sunsetdark,
)

fig.update_layout(
    title_text="Average concentration of NO2 across all years",
)

fig.show()

Dla NO2, rozkład pomiędzy krajami jest już nieco bardziej równomiarny, ze zdecydowanie wyróżniającym się Iranem, a także Mongolią i Chinami.

City-based analysis¶

Correlation (point-biserial)¶

In [ ]:
cities_corr_df = data.copy()
cities_corr_df["is_India"] = cities_corr_df["WHO Country Name"].apply(
    lambda x: 1 if x == "India" else 0
)
cities_corr_df["is_China"] = cities_corr_df["WHO Country Name"].apply(
    lambda x: 1 if x == "China" else 0
)
cities_corr_df
Out[ ]:
WHO Region ISO3 WHO Country Name City or Locality Measurement Year PM2.5 PM10 NO2 PM2.5 temporal coverage PM10 temporal coverage NO2 temporal coverage Reference Number and type of monitoring stations Total number of stations is_India is_China
Loading... (need help?)

Z racji tego, że spodziewano się, iż miasta chińskie oraz indyjskie będą wyróżniały się większym zanieczyzczeniem, uwtorzono zmienne kategoryczne mówiące o przynależności do nich.

In [ ]:
for i in ["is_India", "is_China"]:
    print(f"\n{i}:")
    for c in corr_cols:
        temp_df = cities_corr_df.dropna(subset=c)
        print(f"{i} vs {c}: {stats.pointbiserialr(temp_df[i], temp_df[c])}")
is_India:
is_India vs Measurement Year: SignificanceResult(statistic=-0.033722921784423314, pvalue=1.4305626692267871e-09)
is_India vs PM2.5: SignificanceResult(statistic=0.16586006739676543, pvalue=2.8436607862728626e-93)
is_India vs PM10: SignificanceResult(statistic=0.5670443433607244, pvalue=0.0)
is_India vs NO2: SignificanceResult(statistic=0.03616984654129489, pvalue=7.0204822265107e-08)
is_India vs PM2.5 temporal coverage: SignificanceResult(statistic=-0.3074306764651945, pvalue=4.829027810058524e-159)
is_India vs PM10 temporal coverage: SignificanceResult(statistic=-0.15350792280808412, pvalue=9.846775449244917e-30)
is_India vs NO2 temporal coverage: SignificanceResult(statistic=-0.33641060595561156, pvalue=0.0)
is_India vs Total number of stations: SignificanceResult(statistic=0.08462260994946087, pvalue=2.148838202164032e-15)

is_China:
is_China vs Measurement Year: SignificanceResult(statistic=0.22373939837460843, pvalue=0.0)
is_China vs PM2.5: SignificanceResult(statistic=0.6725632858954628, pvalue=0.0)
is_China vs PM10: SignificanceResult(statistic=0.06486150067896079, pvalue=3.98289492759937e-21)
is_China vs NO2: SignificanceResult(statistic=0.14946638223898057, pvalue=4.4392681451651326e-111)
is_China vs PM2.5 temporal coverage: SignificanceResult(statistic=0.1459973168049667, pvalue=5.979949758409239e-36)
is_China vs PM10 temporal coverage: SignificanceResult(statistic=nan, pvalue=nan)
is_China vs NO2 temporal coverage: SignificanceResult(statistic=nan, pvalue=nan)
is_China vs Total number of stations: SignificanceResult(statistic=nan, pvalue=nan)
c:\Users\adamp\anaconda3\envs\ml_moje\lib\site-packages\scipy\stats\_stats_py.py:4424: ConstantInputWarning:

An input array is constant; the correlation coefficient is not defined.

c:\Users\adamp\anaconda3\envs\ml_moje\lib\site-packages\scipy\stats\_stats_py.py:4424: ConstantInputWarning:

An input array is constant; the correlation coefficient is not defined.

c:\Users\adamp\anaconda3\envs\ml_moje\lib\site-packages\scipy\stats\_stats_py.py:4424: ConstantInputWarning:

An input array is constant; the correlation coefficient is not defined.

Wartości współczynników korelacji punktowo-dwuseryjnej pokazały, że gdy dane miasto należy do Indii, istnieje duża szansa, że jest bardzo zanieczyszczone cząsteczkami PM10. W przypadku, gdy miasto należy do Chin, będzie ono najprawdopodobniej wysokie stężenie PM2.5.

Most and least polluted cities (across all years)¶

PM2.5¶
In [ ]:
grouped_city_pm25 = data.groupby(["WHO Country Name", "City or Locality"])[
    "PM2.5"
].mean()
grouped_city_pm25 = pd.DataFrame(grouped_city_pm25).sort_values(by="PM2.5").dropna()
grouped_city_pm25.index[:10].to_list()
Out[ ]:
[('Portugal', 'Porto Santo'),
 ('Norway', 'Todalen'),
 ('Sweden', 'Bredkalen'),
 ('Australia', 'Emu River'),
 ('Australia', 'St Helens'),
 ('Iceland', 'Suðurnesjabær'),
 ('Iceland', 'Husavik'),
 ('Australia', 'Fingal'),
 ('Norway', 'Birkeland'),
 ('Canada', 'Steeper')]
In [ ]:
grouped_city_pm25.index[-10:].to_list()
Out[ ]:
[('China', 'Hetian Shi'),
 ('Pakistan', 'Rawalpindi'),
 ('China', 'Hotan'),
 ('China', 'Kashi Shi'),
 ('India', 'Delhi'),
 ('India', 'Hapur'),
 ('India', 'Noida'),
 ('India', 'Agra'),
 ('Afghanistan', 'Kabul'),
 ('Cameroon', 'Bamenda')]
In [ ]:
grouped_city_pm25 = data.groupby(["WHO Country Name", "City or Locality"])[
    "PM2.5"
].median()
grouped_city_pm25 = pd.DataFrame(grouped_city_pm25).sort_values(by="PM2.5").dropna()
grouped_city_pm25.index[:10].to_list()
Out[ ]:
[('Portugal', 'Porto Santo'),
 ('Norway', 'Todalen'),
 ('Sweden', 'Bredkalen'),
 ('Australia', 'Emu River'),
 ('Australia', 'St Helens'),
 ('Iceland', 'Suðurnesjabær'),
 ('Australia', 'Fingal'),
 ('Canada', 'Steeper'),
 ('Canada', 'Labrador City'),
 ('Iceland', 'Husavik')]
In [ ]:
grouped_city_pm25.index[-10:].to_list()
Out[ ]:
[('China', 'Hetian'),
 ('China', 'Hetian Shi'),
 ('Pakistan', 'Rawalpindi'),
 ('China', 'Hotan'),
 ('India', 'Agra'),
 ('India', 'Delhi'),
 ('India', 'Hapur'),
 ('India', 'Noida'),
 ('Afghanistan', 'Kabul'),
 ('Cameroon', 'Bamenda')]

Podobnie jak w przypadku analizy na podstawie krajów, tak i tutaj przedstawiono najbardziej i najmniej zanieczyszczone miasta badając i średnią i medianę zanieczyszczeń. Widać, że najczystsze pod względem PM2.5 są miasta na północy Europy oraz wybrane lokalizacje w Australii i Kanadzie. Najbardziej zanieczyszczone są miasta chińskie, indyjskie oraz Bamenda w Kamerunie, stolica Afganistanu - Kabul oraz pakistańskie Rawalpindi.

PM10¶
In [ ]:
grouped_city_pm10 = data.groupby(["WHO Country Name", "City or Locality"])[
    "PM10"
].mean()
grouped_city_pm10 = pd.DataFrame(grouped_city_pm10).sort_values(by="PM10").dropna()
grouped_city_pm10.index[:10].to_list()
Out[ ]:
[('Switzerland', 'Jungfraujoch'),
 ('Iceland', 'Suðurnesjabær'),
 ('Norway', 'Todalen'),
 ('Finland', 'Muonio'),
 ('Sweden', 'Bredkalen'),
 ('Germany', 'Garmisch-Partenkirchen'),
 ('Norway', 'Hurdal Municipality'),
 ('Bahamas', 'Nassau'),
 ('Norway', 'Birkeland'),
 ('Sweden', 'Malmberget')]
In [ ]:
grouped_city_pm10.index[-10:].to_list()
Out[ ]:
[('Egypt', 'Greater Cairo'),
 ('Bahrain', 'Ras Hayan'),
 ('Iran (Islamic Republic of)', 'Boshehr'),
 ('India', 'Jharia'),
 ('Bahrain', 'Hamad Town'),
 ('Saudi Arabia', 'Al Jubail'),
 ('Iraq', 'Basra'),
 ('Pakistan', 'Rawalpindi'),
 ('Iran (Islamic Republic of)', 'Zabol'),
 ('Pakistan', 'Peshawar')]
In [ ]:
grouped_city_pm10 = data.groupby(["WHO Country Name", "City or Locality"])[
    "PM10"
].median()
grouped_city_pm10 = pd.DataFrame(grouped_city_pm10).sort_values(by="PM10").dropna()
grouped_city_pm10.index[:10].to_list()
Out[ ]:
[('Switzerland', 'Jungfraujoch'),
 ('Iceland', 'Suðurnesjabær'),
 ('Norway', 'Todalen'),
 ('Finland', 'Muonio'),
 ('Sweden', 'Bredkalen'),
 ('Germany', 'Garmisch-Partenkirchen'),
 ('Norway', 'Hurdal Municipality'),
 ('Bahamas', 'Nassau'),
 ('Sweden', 'Malmberget'),
 ('Italy', 'Alpe Devero')]
In [ ]:
grouped_city_pm10.index[-10:].to_list()
Out[ ]:
[('Bahrain', 'Ras Hayan'),
 ('Iran (Islamic Republic of)', 'Boshehr'),
 ('Egypt', 'Greater Cairo'),
 ('India', 'Jharia'),
 ('Bahrain', 'Hamad Town'),
 ('Saudi Arabia', 'Al Jubail'),
 ('Iraq', 'Basra'),
 ('Pakistan', 'Rawalpindi'),
 ('Iran (Islamic Republic of)', 'Zabol'),
 ('Pakistan', 'Peshawar')]

Dla PM10 wyniki prezentują się podobnie, dodając do "czystych miast" Jungfraujoch w Szwajcarii, Garmisch-Partenkirchen w Niemczech oraz Nassau na Bahamach. Najbardziej zanieczyszczone są miasta leżące całkiem blisko siebie w Azji Południowo-Wschodniej oraz Egipcie.

NO2¶
In [ ]:
grouped_city_no2 = data.groupby(["WHO Country Name", "City or Locality"])["NO2"].mean()
grouped_city_no2 = pd.DataFrame(grouped_city_no2).sort_values(by="NO2").dropna()
grouped_city_no2.index[:10].to_list()
Out[ ]:
[('Canada', 'Aylesford'),
 ('Canada', 'Southampton'),
 ('Switzerland', 'Jungfraujoch'),
 ('Norway', 'Hemnes'),
 ('Sweden', 'Bredkalen'),
 ('Norway', 'Todalen'),
 ('Austria', 'Rauris'),
 ('Italy', "Sant'eufemia A Maiella"),
 ('Italy', 'Seulo'),
 ('Finland', 'Muonio')]
In [ ]:
grouped_city_no2.index[-10:].to_list()
Out[ ]:
[('Mexico', 'Zona Metropolitana Del Valle De Toluca'),
 ('Brazil', 'Belo Horizonte'),
 ('Kuwait', 'Ali Subah Al-Salem'),
 ('Iran (Islamic Republic of)', 'Tehran'),
 ('Iran (Islamic Republic of)', 'Varamin'),
 ('Iran (Islamic Republic of)', 'Orumye'),
 ('Iran (Islamic Republic of)', 'Pakdasht'),
 ('Iran (Islamic Republic of)', 'Hamedan'),
 ('Iran (Islamic Republic of)', 'Yasoj'),
 ('Iran (Islamic Republic of)', 'Arak')]
In [ ]:
grouped_city_no2 = data.groupby(["WHO Country Name", "City or Locality"])[
    "NO2"
].median()
grouped_city_no2 = pd.DataFrame(grouped_city_no2).sort_values(by="NO2").dropna()
grouped_city_no2.index[:10].to_list()
Out[ ]:
[('Canada', 'Aylesford'),
 ('Canada', 'Southampton'),
 ('Switzerland', 'Jungfraujoch'),
 ('Norway', 'Hemnes'),
 ('Sweden', 'Bredkalen'),
 ('Norway', 'Todalen'),
 ('Austria', 'Rauris'),
 ('Italy', "Sant'eufemia A Maiella"),
 ('Italy', 'Seulo'),
 ('Finland', 'Muonio')]
In [ ]:
grouped_city_no2.index[-10:].to_list()
Out[ ]:
[('Mexico', 'Zona Metropolitana Del Valle De Toluca'),
 ('Brazil', 'Belo Horizonte'),
 ('Kuwait', 'Ali Subah Al-Salem'),
 ('Iran (Islamic Republic of)', 'Varamin'),
 ('Iran (Islamic Republic of)', 'Tehran'),
 ('Iran (Islamic Republic of)', 'Orumye'),
 ('Iran (Islamic Republic of)', 'Pakdasht'),
 ('Iran (Islamic Republic of)', 'Hamedan'),
 ('Iran (Islamic Republic of)', 'Yasoj'),
 ('Iran (Islamic Republic of)', 'Arak')]

W przypadku NO2 także przodują miasta północy, dodając do nich austriackie Rauris, a także włoskie Sant'eufemia A Maniella i Seulo. Najbardziej zanieczyszczone są praktycznie wszystkie miasta Iranu oraz Ali Subah Al-Salem w Kuwejcie, Belo Horizonte w Brazylii oraz Zona Metropolitana Del Valle De Toluca w Meksyku.

Biggest progress and regress in air polution¶

PM2.5¶
In [ ]:
grouped = data.groupby(
    ["WHO Country Name", "City or Locality", "Measurement Year"]
).mean(numeric_only=True)
grouped_pm25 = grouped.dropna(subset="PM2.5")

earliest = grouped_pm25.groupby(level=(0, 1))["PM2.5"].first()
latest = grouped_pm25.groupby(level=(0, 1))["PM2.5"].last()

difference = earliest - latest  # + --> progress; - --> regres

difference
Out[ ]:
PM2.5
WHO Country Name City or Locality
Loading... (need help?)
In [ ]:
difference.sort_values()[:10]
Out[ ]:
PM2.5
WHO Country Name City or Locality
Loading... (need help?)
In [ ]:
difference.sort_values()[-10:]
Out[ ]:
PM2.5
WHO Country Name City or Locality
Loading... (need help?)

Przechodząc do analizy progresu i regresu, widać, że w przypadku PM2.5 po obu stronach barykady znajdują się miasta chińskie, w tym sama stolica. Największy postęp dokonał się w irańskim Zabol, a zwiększenie zanieczyszczenia w pakistańskim Lahore.

PM10¶
In [ ]:
grouped_pm10 = grouped.dropna(subset="PM10")

earliest = grouped_pm10.groupby(level=(0, 1))["PM10"].first()
latest = grouped_pm10.groupby(level=(0, 1))["PM10"].last()

difference = earliest - latest  # + --> progress; - --> regres

difference
Out[ ]:
PM10
WHO Country Name City or Locality
Loading... (need help?)
In [ ]:
difference.sort_values()[:10]
Out[ ]:
PM10
WHO Country Name City or Locality
Loading... (need help?)
In [ ]:
difference.sort_values()[-10:]
Out[ ]:
PM10
WHO Country Name City or Locality
Loading... (need help?)

W przypadku PM10 ogromny regres zaszedł w Pasakha w Bhutanie, w Ma'ameer w Bahrajnie i Gorakpurze w Indiach, a progres: w stolicu Iraku, w Qazvin w Iranie oraz w Nanded w Indiach.

NO2¶
In [ ]:
grouped_no2 = grouped.dropna(subset="NO2")

earliest = grouped_no2.groupby(level=(0, 1))["NO2"].first()
latest = grouped_no2.groupby(level=(0, 1))["NO2"].last()

difference = earliest - latest  # + --> progress; - --> regres

difference
Out[ ]:
NO2
WHO Country Name City or Locality
Loading... (need help?)
In [ ]:
difference.sort_values()[:10]
Out[ ]:
NO2
WHO Country Name City or Locality
Loading... (need help?)
In [ ]:
difference.sort_values()[-10:]
Out[ ]:
NO2
WHO Country Name City or Locality
Loading... (need help?)

Dla wartości NO2 zaszły mniejsze zmiany niż dla PM10 oraz PM2.5, ale największy regres zanotowano głównie w miastach indyjskich, a progres w Varamin w Iranie, Kagoszimie w Japoni oraz stolicy Iranu.

Stations analysis¶

In [ ]:
summed_by_cities_stations_df = (
    data.groupby(["WHO Country Name", "City or Locality"])["Total number of stations"]
    .sum()
    .reset_index()
)

summed_by_cities_stations_df.sort_values("Total number of stations", ascending=False)[
    :10
]
Out[ ]:
WHO Country Name City or Locality Total number of stations
Loading... (need help?)

Tutaj zbadano, w których miastach na świecie znajduje się najwięcej stacji pomiarowych wedle dostępnych danych. Jak widać najwięcej jest ich w Kalkucie, Abu Dhabi oraz Pradze.

Air pollution map¶

In [ ]:
# from geopy.geocoders import Nominatim
# from geopy.exc import GeocoderTimedOut

# geolocator = Nominatim(user_agent="EDA")

# data_lon_lat = data.copy()

# cities_list = data_lon_lat["City or Locality"].unique()
# cities_df = pd.DataFrame({"City": cities_list, "Longitude": None, "Latitude": None})


# def get_coordinates(city):
#     try:
#         location = geolocator.geocode(city)
#         if location:
#             return location.longitude, location.latitude
#         else:
#             return None, None
#     except GeocoderTimedOut as e:
#         print("Error: geocode failed on input %s with message %s" % (city, e.message))


# cities_df[["Longitude", "Latitude"]] = cities_df["City"].apply(
#     lambda x: pd.Series(get_coordinates(x))
# )
# cities_df

W celu wykonania mapy miast i zanieczyszczeń próbowano zmapować miasta na współrzędne geograficzne, co ostatecznie się nie udało ze względu na zbyt dużą liczbę zapytań (zbyt dużą liczbę miast w analizowanym zbiorze danych).

W związku z tym wspomożono się plikiem dostępnym pod tym adresem: https://simplemaps.com/data/world-cities

Get cities coords¶
In [ ]:
cities_gdf = gpd.read_file("../data/worldcities.csv")

data_lat_lng = data.merge(cities_gdf, left_on="City or Locality", right_on="city")
print(len(data_lat_lng["City or Locality"].unique()))
print(len(data["City or Locality"].unique()))
data_lat_lng
2535
6874
Out[ ]:
WHO Region ISO3 WHO Country Name City or Locality Measurement Year PM2.5 PM10 NO2 PM2.5 temporal coverage PM10 temporal coverage NO2 temporal coverage Reference Number and type of monitoring stations Total number of stations city city_ascii lat lng country iso2 iso3 admin_name capital population id geometry
Loading... (need help?)

W pliku tym znajdują się wszystkie stolice państw oraz stanów/województw (prawie połowa z analizowanych), zatem uznano ten zbiór za reprezentatywny.

PM2.5¶
In [ ]:
grouped_city_pm25 = data_lat_lng.groupby(
    ["WHO Country Name", "City or Locality", "lat", "lng"]
)["PM2.5"].mean()
df_grouped_city_pm25 = grouped_city_pm25.reset_index().dropna()

fig = px.scatter_geo(
    df_grouped_city_pm25,
    lat=df_grouped_city_pm25["lat"],
    lon=df_grouped_city_pm25["lng"],
    color="PM2.5",
    color_continuous_scale=px.colors.sequential.Sunsetdark,
    hover_name="City or Locality",
)

fig.update_layout(
    title_text="Average concentration of PM2.5 across all years",
)

fig.show()

Utworzona mapa zanieczyszczeń poszczególnych miast pokazuje niektóre miasta Afryki oraz całą Azję Południową jako najbardziej zanieczyszczone, a Amerykę Północną, Australię oraz Europę Północną jako te z najlepszym powietrzem pod względem PM2.5.

PM10¶
In [ ]:
grouped_city_pm10 = data_lat_lng.groupby(
    ["WHO Country Name", "City or Locality", "lat", "lng"]
)["PM10"].mean()
df_grouped_city_pm10 = grouped_city_pm10.reset_index().dropna()

fig = px.scatter_geo(
    df_grouped_city_pm10,
    lat=df_grouped_city_pm10["lat"],
    lon=df_grouped_city_pm10["lng"],
    color="PM10",
    color_continuous_scale=px.colors.sequential.Sunsetdark,
    hover_name="City or Locality",
)

fig.update_layout(
    title_text="Average concentration of PM10 across all years",
)

fig.show()

Mapa PM10 prezentuje się analogicznie do mapy z PM2.5. Tutaj jedynie ze względu na dostępność danych o PM10 oraz o miastach, to głównie Indie i Pakistan są wskazane jako najbardziej zanieczyszczone.

NO2¶
In [ ]:
grouped_city_no2 = data_lat_lng.groupby(
    ["WHO Country Name", "City or Locality", "lat", "lng"]
)["NO2"].mean()
df_grouped_city_no2 = grouped_city_no2.reset_index().dropna()

fig = px.scatter_geo(
    df_grouped_city_no2,
    lat=df_grouped_city_no2["lat"],
    lon=df_grouped_city_no2["lng"],
    color="NO2",
    color_continuous_scale=px.colors.sequential.Sunsetdark,
    hover_name="City or Locality",
)

fig.update_layout(
    title_text="Average concentration of NO2 across all years",
)

fig.show()

W przypadku NO2 nie można określić konkretnego regionu jako najbardziej zanieczyszczonego, może oprócz Teheranu w Iranie. Tutaj niestety wartości NO2 są wszędzie podobnie wysokie.

Na pytania 1-4, a także częściowo kolejne, wypowiedziano się w trakcie analizy. Poniżej odpowiedzi na pytania 5-7.

  1. Na podstawie przeprowadzonej analizy widać wyraźnie dodatnią zależność pomiędzy wartościami PM2.5 oraz PM10, co jest zgodne z codzienymi obserwacjami - zanieczyszczenia ta powstają w wyniku emisji spalin samochodowych oraz paleniu w piecach, zatem oba są produkowane w mniej więcej podobny sposób. Ponadto można wskazać rejony Bliskiego Wschodu oraz Azji Południowo-Wschodniej jako te najbardziej zanieczyszczone. Co ważne, często (choć nie zawsze) są to regiony biedne, kraje rozwijające się, co można dość łatwo połączyć. W takich miejscach panuje zdecydowanie mniejsza świadomość ekologiczna, jest mniejszy współczynnik ludzi wykształconych, a żyjący tam ludzie mają często po prostu inne priorytety. W związku z tym regiony te od lat są zanieczyszczone, co koreluje się z niskim poziomem życia.
  2. Ten rodzaj danych może mieć jedynie bardzo ogólne znaczenie jako dane wejściowe do modelu ML. Dane te są zagregowane do roku, co jest dość dużym stopniem generalizacji, ponadto występuje tu wiele braków i można by przewidzieć jedynie ogólne trendy, tzn. które miejsca ogólnie staną się jeszcze bardziej, lub też nieco mniej zanieczyszczone. Należałoby dodać choćby dane miesięczne, aby np. móc predykować w jakich okresach roku należy nosić maseczki antysmogowe w danych miastach, czy też wprowadzać inne obostrzenia. Niewątpliwie jednak dane te mogą pomoć wytworzyć pewien ogólny model zanieczyszczeń.
  3. Przedstawiona eksploracyjna analiza danych pokazuje jak wiele można uzyskać z praktycznie surowych danych. Jesteśmy w stanie zobaczyć trendy zanieczyszczeń, ocenić najbardziej oraz najmniej zanieczyszczone miejsca, a także określić relacje między zmiennymi w oparciu o różne poziomy generalizacji. Dowiedziono, że Europa Północna słynie z czystego powietrza, podobnie jak Kanada oraz Australia i Nowa Zelandia. Ponadto przedstawiono rejon Iranu oraz krajów sąsiadujących (także "przez morze") jako te najbardziej skażone dwutlenkiem azotu. Również ukazano, że Indie oraz Chiny są świadome stopnia swojego zanieczyszczenia cząsteczkami PM2.5 oraz PM10 - to tam znajduje się najwięcej czujników jednocześnie wskazujących najwyższe zanieczyszczenie. Widać także, że oba te kraje starają się polepszyć sytuację w niektórych miastach, co przedstawia różnica w zanieczyszczeniach na korzyść ostatnich lat dla wybranych chińskich i indyjskich miast. Tak przeprowadzone EDA, pomaga ocenić, które zmienne uwzględnić w modelu, oraz który stopień generalizacji przyjąć. Pomaga także w ocenie danych wyjściowych z modelu - po porównaniu ich z tego typu analizą danych historycznych możemy ogólnie ocenić, czy model działa dobrze, czy też nie.